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Abstract 

A key step in the application of formal automated design techniques to structures un- 
der transient loading is the calculation of sensitivities of response quantities to the design 
parameters. This paper considers structures with general forms of damping acted on by 
general transient loading and addresses issues of computational errors and computational 
efficiency. The equations of motion are reduced using the traditional basis of vibration 
modes and then integrated using a highly accurate, explicit integration technique. A crit- 
ical point constraint formulation is used to place constraints on the magnitude of each 
response quantity as a function of time. Three different techniques for calculating sensi- 
tivities of the critical point constraints are presented. The first two techniques are based 
on the straightforward application of the forward and central difference operators, respec- 
tively. The third technique is based on explicit differentiation of the equations of motion. 
Condition errors, finite difference truncation errors, and modal convergence errors for the 
three techniques are compared by applying them to a simple five-span-beam problem. 
Sensitivity results are presented for two different transient loading conditions and for both 
damped and undamped cases. 

* Aerospace Engineer, Structures and Dynamics Division, 
f Professor, Department of Aerospace and Ocean Engineering. 
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Nomenclature 


C = system damping matrix 

C = reduced damping matrix 

f = vector of external forces 

? = vector of reduced system forces 

<7,- = *th constraint 

K = system stiffness matrix 

K = reduced stiffness matrix 

M = system mass matrix 

M = reduced mass matrix 

m = number of equations in the reduced system 
N,y = submatrix i j in matrix series expansion method 
n = number of degrees-of-freedom in the finite element model 
q = vector of reduced system coordinates 
t = time 

t c = critical time in critical point constraint 

u = vector of displacements 

u i = displacement at the one-third point in the *th span of the five-span-beam example 
u a iiow — maximum allowable value of a displacement 

Xi — the »th design variable 

x = vector of design variables 

W ij = submatrix i j in matrix series expansion method 
$ = matrix with each column being a basis vector 

Introduction 

A key step in most automated structural design approaches is the calculation of deriva- 
tives of the structural behavior with respect to the design variables (i.e. sensitivities). This 
problem has been addressed in references describing overall automated design procedures 
and in references which focus on procedures for sensitivity calculations. Reference [l] 
provides a comprehensive review of work on calculating sensitivities in structural systems. 
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Most of this work has addressed static structural response and structural behavior 
governed by eigenvalue problems, but a few researchers have considered the problem of 
sensitivity calculations and optimization in transient structural dynamics problems. Ref- 
erence [2] considers the minimum mass design of structures under displacement and stress 
constraints but assumes a simplified form of external, transient loading in order to elimi- 
nate the time dependence. Reference [3] considers only harmonic loading but concentrates 
on the difficult problem of the disjoint design space that often occurs in the automated 
design problem under dynamic loading. 

Another class of work has considered the application of general sensitivity theory to 
the problem of dynamic mechanical systems under arbitrary loading. Reference [l] sum- 
marizes this work and describes three basic approaches which have been employed. In the 
first method, called the direct method, the transient equations are directly differentiated 
and solved. In the second method, called the adjoint method, derivatives of a constraint 
functional are obtained with the aid of an adjoint variable. In the third method, called 
the Green’s function method, the derivatives are obtained in terms of the Green’s function 
of the equations of motion. Although the results from all three methods are theoretically 
identical, their relative computational efficiency depends on the relative numbers of design 
variables, degrees of freedom, and constraints. Haug and his co-workers have made con- 
siderable progress in developing the adjoint variable method and applying it to problems 
in machine design and elastic structural design (e.g. refs. [4], [5]). The machine design 
problems often have the additional complexity of nonlinear equations of motion. However 
most of these example problems have involved few degrees-of-freedom, design variables, 
or constraints. There has been little application of these transient sensitivity calculation 
methods to large-order structural design problems where computational issues become 
especially important. 

The objective of this paper is to address these computational issues in calculating sen- 
sitivities in transient structural dynamics problems. In addition to the analytically-based 
sensitivity formulations described above, straightforward finite difference procedures are 
also strong candidates for use in transient structural problems. Both a single-step, forward 
difference procedure and a two-step, central difference procedure are considered herein. A 
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version of the direct analytical method described above is also considered. This method 
uses the directly differentiated transient equations but calculates the derivatives of the co- 
efficients in the transient equations using a forward difference operator. As a consequence, 
the method is termed semi-analytical. In addition to the selection of sensitivity calculation 
method, several other computational issues are addressed herein. The first issue concerns 
numerical solution of the transient equations and the sensitivity equations in the semi- 
analytical method. Any numerical error associated with this solution process could have a 
significant effect on the accuracy of derivatives. A second issue is the particular constraint 
formulation used to represent the time-dependent constraints. Herein, a formulation is 
used where constraints are enforced only at a few, selected, critical points in time. A pro- 
cedure for selecting these critical times is developed. A third issue concerns the selection of 
finite difference step size on the accuracy of derivatives. This is a key issue whenever finite 
difference operators are used and is important in all three sensitivity methods considered 
herein. A fourth issue is the convergence of the modal method, traditionally used to re- 
duce the computational cost in large-order, linear dynamics problems, as a function of the 
number of modes used. The particular concern here is the convergence of the sensitivities. 
Finally, the last issue is the computational cost of the sensitivity calculation methods. This 
is discussed in general, relative terms rather than with specific performance measures. 

The investigation of the above issues is carried out by applying the two finite difference 
methods and the semi-analytical method to a five-span-beam example. In some cases, 
concentrated, point masses and dampers are added to the beam. Two different loading 
conditions are also considered. 

Equations of Motion and Solution 

Governing Equations 

The equations of motion for a damped structural system can be written as 

Mu + Cu + Ku = f (t) (1) 
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which is a set of n, coupled differential equations. These equations are usually of high 
order and are typically stiff. Therefore integration techniques used to solve them are 
typically unconditionally stable and produce an accurate solution only for the low frequency 
components. The well-known Newmark method (see for example ref. [6]) is an example of 
such an integration technique. 

An alternate approach to the direct integration of eq. (1) is to express the solution in 
terms of m basis functions 

u = $q (2) 

where m is usually much less than n. Then a reduced set of equations can be written 

I&q + Cq + Kq = ? (t) (3) 

where M = 9 T M9, C = * T C$, K = S r KS, and f = $ T f. 

When the undamped natural vibration modes of the system are used as basis functions 
and C is a linear combination of M and K , equation (3) represents a set of uncoupled 
differential equations which can be solved independently. For certain forms of external 
loading, such as a piecewise linear approximation, an exact, explicit solution equation is 
available. This approach is described in [7] and is used in NASTRAN [8]. 

When vectors which are not the vibration mode shapes are used as basis functions, or 
a more general form of damping is modeled, the dynamic equations are coupled and must 
be solved simultaneously. In contrast to eqs. (1), however, eqs. (3) are low order, not stiff, 
and the primary concern is accurately integrating every equation in the system. Therefore 
an integration method which minimizes truncation errors in the solution is preferred even 
if it is only conditionally stable. Accuracy is especially important in sensitivity analyses 
because errors in the solution process are usually magnified in the calculation of derivatives. 
Eliminating the integration error due to time step selection can significantly improve the 
accuracy of sensitivities. 
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Solution Method 

An approach that allows the use of moderately large time steps and effectively elimi- 
nates truncation error is called the matrix series expansion method in [9] and the transfer 
matrix method in [10] and [11]. This method expands the solution in a Taylor series where 
the number of terms determines the accuracy of the approximation. Using this series an 
expression can be written for the solution at time t + At in terms of the solution and load 
at time t. 



It has been assumed here that the load f is approximated as a piecewise linear function 
of time and therefore the second and higher order derivatives equal zero. The values of 
the coefficients W,y and N»y depend on the number of terms taken in the series. Since 
the number of terms in the series has only a very small effect on the computational cost 
of the method, 50 terms were used in this study to effectively eliminate truncation error. 
Another consideration is that the matrix series expansion method is only conditionally 
stable and therefore a critical value of At exists. Numerical experimentation has shown 
that At must be smaller than the shortest period of the system to insure numerical stability. 
Unconditionally stable, low-order methods such as the Newmark method would require 
much smaller time steps for comparable accuracy. 

Critical Point Constraint 


Constraint Formulation 

The general form of the constraint equation is 


0 ,(x,t)<O (5) 

where x is a vector of design variables and t is time. An effective approach for insuring 
that this constraint is satisfied for all values of t is the “critical point constraint” approach 
described in ref. [ 12 ] pages 168-169. In this approach a set of peak values of the function 
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g (denoted critical points) is selected. An obvious point to include is the time with the 
“worst* value of g. However, if only this point is included, an optimization process modi- 
fying a structure based on this information might unknowingly produce a design where the 
constraint is violated at another time point. To guard against this possibility, a number 
of important peaks are selected. References [13] and [14] consider in detail the efficient 
location of critical points in large-scale structures problems with many constraints. The 
present paper presents a method for selecting the most important peaks as critical points 

In the work reported herein, constraints are placed on selected displacements such 

that 

Ju t '(x,l)| < Mallow (6) 

where are the displacements at specific points in the structure and u a u ow is the allowable 
value of the displacement. The critical values of this constraint occur at points in time 
where u t - has the largest magnitude. These are identified by examining every value of u 
along the time history. In the implementation here, each constraint is assumed to have a 
specified number of critical points; five critical points for each u,- are selected. Values of u 
where du/dt = 0 or values of u at the end points of the time interval are local maxima of 
g and are termed candidate critical points. 

Selection of Critical Points 

The procedure for selecting the critical points from these candidates can best be 
explained by referring to an example displacement time history shown in figure 1. The 
critical points are labeled with circled numbers and a few of the many candidate critical 
points ^e labeled with circled letters. The selection criteria applied to every candidate 
critical point will be explained by considering these few candidate points. Candidate 
critical points a and c were discarded because the absolute values of the displacements 
at these points were smaller than those at the five other critical points. The criterion 
for discarding candidate points b, d, and e is slightly more complicated. From figure 1 
it can be seen that all three candidate points have larger displacement magnitudes than 
that of critical point 1, for example. However, candidate points b, d, and e are all part of 
“major” peaks where a critical point is selected. A second criterion applied to the selection 
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process is a requirement that only one critical point from each major peak be selected. This 
insures that the critical points represent the total dynamic response rather than just the 
high frequency undulations on, at worst, a single major peak. 

A major peak is identified with the following procedure. Whenever a critical point is 
selected after comparing its magnitude with that at other critical points, a special screening 
process is activated. This screening process tests the displacement at every subsequent time 
point to determine if it differs from that at this last selected critical point by at least a 
specified percentage (25% for the studies reported herein ). If so, all subsequent time 
points are no longer considered part of the current major peak. Any candidate critical 
points identified while this special screening process is in effect are compared only against 
the last selected critical point. 

An example is the major peak in figure 1 which contains points d and 4. In the 
selection process, d is initially selected as a critical point and the screening process is 
activated. The three points where du/dt = 0 between point d and point 4 are recognized 
to be part of the same major peak as d, but since the magnitude of the displacements at 
these points is smaller than at point d, they are discarded. Point 4 is also part of the same 
major peak as d but since the displacement magnitude there is larger than at point d, it 
replaces d as a critical point. Before the next candidate critical point is considered, the 
displacement has changed from that at point 4 by more than 25% and so is considered to 
be on a new major peak. 

Derivatives of Critical Point Constraints 

Once the critical points have been identified for the nominal design, these can used 
in calculating sensitivities. Reference [12] demonstrates that the change in time location 
of critical points can be neglected in calculating derivatives of peak values with respect to 
design variables by examining the expression for the total derivative of g with respect to 
a design variable x. Considering a constraint g(x, t) at a critical time t c 



The last term in (7) is always zero because at interior critical points dg/dt = 0 and on 
the boundary dt e /dx = 0. Accordingly, the sensitivity calculations need to be performed 
only at the specific times where the critical points have been identified. This can result in 
a considerable savings in computational time, especially when there are many constraints, 
many time points, or many basis functions used to represent the response. The details of 
each sensitivity calculation method are discussed below. 

Methods For Calculating Sensitivities 


Finite Difference Methods 

Both the forward difference and central difference methods have been used in this 
study to calculate sensitivities. The well known forward difference approximation to du/dx, 


Au _ u(x + Ax) — u(x) 
Ax Ax 


and central difference approximation 


( 8 ) 


Au _ u(x + Ax) — u(x — Ax) 

Ax 2Ax 

are used. The truncation error for the forward difference approximation is 

«r(Ax) = -^^(x + fAx) 0<f <1 

and is 


(9) 


( 10 ) 


«t(Ax) = -^-^(x + ;Ax) 0<c<l (11) 

for the central difference approximation. 

In applying equations (8) and (9), the selection of difference step size Ax is a concern. 
Selection of a large step size results in errors in the derivative due to truncation of the 
operator (eqs. 10 and 11). Selection of a small step size can lead to errors in the derivative 
due to the limited floating point precision of the computer or algorithmic inaccuracies in 
calculating u (condition errors). Frequently, no acceptable value exists for Ax to produce 
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an accurate value of du/dx considering the conflicting requirements of minimizing trun- 
cation and condition errors. Because the truncation error associated with the operator of 
equation 9 is typically less that that of equation 8, it is possible to use a larger finite dif- 
ference step size. The larger Ax reduces the condition error from the function evaluations 
and results in a more accurate value of du/dx. However, the necessity of two function 
evaluations needed for equation 9 makes the procedure more computationally costly. 

For the studies in this paper the natural vibration mode shapes are used as basis 
functions to represent the transient response. In calculating the response of the perturbed 
design in equation (8) and the two perturbed designs in equation (9) some computational 
savings are possible relative to the computations for the initial design. 

If the mode shapes for the initial design are used to represent the perturbed design, 
the cost of re-solving the eigenvalue problem is eliminated. However, the reduced set of 
equations for the perturbed system must still be formed and M , C , R are now full. This 
coupled system is then solved using the matrix series expansion method described above. 

If the updated mode shapes for the perturbed design are used in the analysis, many 
eigensolution procedures, such as the subspace iteration used here, can begin with the 
mode shapes from the initial design as approximations. Since the perturbation in the 
design is small, the subspace iteration procedure converges rapidly. However, at least one 
factorization of K is required. For large finite element models this can be the largest part 
of the computational cost. For most of the studies in this paper, the forward difference 
method used the initial mode shapes to represent the perturbed design. Because the central 
difference method was used for reference values of derivatives, updated mode shapes were 
calculated for the two required perturbed designs. In both case, because of the critical point 
constraint formulation, the transformation from modal coordinates to physical coordinates 
(displacements, stresses, etc.) is performed only at the critical times instead of at all time 
points. 
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Semi- Analytical Method 

The direct method for sensitivity calculation is derived by differentiating equation 
(3). This is presented in detail in [12] pages 169-171 and will be summarized here. After 
differentiating (3) with respect to the design variable ,x, the result is 


-dq .dd i -dq df dM dC . dK 

dx dx dx dx dx ^ dx ^ dx ^ 

The first step in forming this equation is the calculation of the derivatives of f , M , C , 

and K (equation 3) with respect to x. This process is simplified considerably by using a 

fixed set of basis functions, $ . For the studies in this paper $ is taken as a set of vibration 

modes of the initial structure. With $ independent of x, differentiation of K yields 



_ 'V uiV _ / rt v 

n • (13) 

with similar expressions for the derivatives of M and C . 

The force f is frequently not a function of the design variables (as in the models 
of this paper) or is an explicit function of x and therefore easy to differentiate. The 
derivatives of M , C , and K are, in general, difficult to calculate and for this reason these 
derivatives are often replaced by finite difference approximations. The combination of 
analytical differentiation of the response equations with finite difference matrix derivatives 
is known as a semi-analytical approach. Semi-analytical methods have become the methods 
of choice in static, structural sensitivity calculations because they allow generality in the 
selection of design variables with moderate computational cost. Herein a semi-analytical 
implementation of the direct method (abbreviated as “the semi-analytical method”) is used 
with the matrix derivatives calculated using the forward difference operator. For several 
important classes of design variables, however, M , C , and K are linear functions of x. 
For example, M and K in a finite element model composed of truss members are linear 
functions of member cross sectional area. In these cases there are no truncation errors 
and large finite difference step sizes can be used to reduce the condition error and produce 
accurate derivatives. 
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After forming the derivatives of ? , M , C , and R the right-hand side of (12) can be 
formed using ^ ^ , and q from the solution of (3). The matrix series expansion method 
insures that accurate values of q , ^ , and q are available for this step. . Equation (12) 
can then be integrated in time to yield dq/dx , di\jdx , and dq/dx. Again, because of the 
critical point constraint formulation, the transformation from these modal quantities to 
physical quantities is performed only at the critical times. 

Computational Cost of the Methods 

All three methods require an identical analysis of the initial design. If updated vibra- 
tion modes are used in the two finite difference methods, this analysis must be repeated 
resulting in a high computational cost for calculating the sensitivities. For large problems, 
the cost of re-solving the vibrational eigenproblem for each perturbed case would become 
excessive. The two methods presented here with reasonable costs for large problems are 
the forward difference method with fixed modes and the semi-analytical method. The 
forward difference method requires M , C , and K for the perturbed design. The semi- 
analytical method, however, requires formation of dM/dx, dC/dx, and dK/dx. The costs 
for these two operations are roughly the same. The forward difference method requires 
the integration of the reduced equations of motion (3) for the perturbed model. The 
semi-analytical method requires integration of (12). The costs for these two operations 
are similar. The semi-analytical method, however, requires the formation of the right- 
hand-side in eq. (12) and there is no similar operation in the forward difference method. 
To form the right-hand-side of (12) three matrix multiplications must be performed at 
each time step. Similar operations are performed in integrating the reduced equations of 
motion (see eq. 4) so the cost of forming the right-hand-side of (12) is of the same order 
as an additional time integration. This cost will not become large unless the order of the 
reduced system, m, or the number of time steps become large. 
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Computer Implementation 

The three methods for calculating sensitivities have been implemented using a slightly 
modified version of the general purpose finite element code, EAL [9]. EAL includes general 
language constructs for controlling execution flow as well as general and specific utilities 
for manipulating data stored as named entities in a database. Each method was imple- 
mented as a set of EAL procedures which are invoked in a manner similar to FORTRAN 
subroutines. 

The transient response module in EAL version 312 solves the uncoupled form of 
equation (3) using the matrix series expansion method. A modification was made to 
allow equation (3) to be fully coupled. In the semi-analytical method, the right-hand 
side, pseudo loading of equation (12) can be easily formed using EAL. However, a slight 
modification to the transient response module was required to permit solution of equation 
(12) with this general form of loading. In addition, a special purpose module was added 
to EAL to perform the task of identifying the critical points on each response function. 

Most of the calculations were performed using a CDC Cyber version of EAL where 
the floating point operations were done with 60 bit precision. In order to demonstrate the 
effects of condition error, a few calculations were performed using a DEC VAX version of 
EAL where most of the floating point operations were done with 32 bit precision. 

Five-span-beam Example 

The three methods for calculating sensitivities axe applied to a five-span, planar beam 
example taken from reference [15] and shown in figure 2. The beam is modeled with three 
beam finite elements per span resulting in 26 unconstrained degrees-of-freedom. As shown 
in figure 2, translational and rotational viscous dampers were also added to the beam. 
These devices are representative of velocity feedback controllers which might be added 
to flexible structures. Cases with and without dampers were considered. The numerical 
values of the damping coefficients from ref. [15] of C\ = .008 sec-lbf/in and c 2 = 1.2 sec-lbf 
were used. A case was also considered where a 1.0 lb mass (approximately 20% of the 
beam’s mass) was added to the beam at the location of the translational damper. The 
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Table 1. Eigenvalues For the Three Five-span-beam Cases 


Mode 

Undamped 

Damped 

Undamped 
With Point Mass 


Frequency, Hz. 

Frequency, Hz. 

Damping Ratio 

Frequency, Hz. 

1 

1.1707 

1.2210 

.0851 

.9401 

2 

1.2991 

1.2926 

.0352 

1.2594 

3 

1.6254 

1.6298 

.0690 

1.5445 

4 

2.0491 

2.0910 

.0590 

1.8005 

5 

2.4628 

2.5497 

.0958 

2.3729 

6 

4.7343 

4.8426 

.0044 

4.2327 

7 

5.0105 

4.9785 

.0413 

4.8858 

8 

5.6472 

5.7703 

.0126 

5.6400 

9 

6.4153 

6.4178 

.0407 

5.9261 

10 

7.1274 

7.2229 

.0193 

6.8762 


eigenvalues for these three cases are shown in table 1. The additional point mass has a 
significant effect on the frequencies while the dampers have little effect. 

Effect of Loading and Damping on Response 

The applied loading for all studies consisted of a point moment of .04405 in-lbs applied 
at the right end of the beam. Two different time functions for this load, a step and a ramp, 
(shown in figure 2) were considered. The step load tends to excite many of the higher modes 
of the beam. This is shown in figure 3a where a selected displacement for the undamped 
beam, u 2 (figure 2), is plotted as a function of time. The participation of the higher 
frequency modes in the response is obvious from the jagged nature of the peaks. However, 
when the load is applied more gradually with the ramp function shown in figure 2, the 
higher modes are excited less. This is shown in figure 3b where displacement u 2 is plotted 
as a function of time. The addition of the two point dampers also serves to reduce the 
participation of the high frequency modes for the beam under the step load. Displacement 
u 2 is plotted as a function of time in figure 3c for this case. All 26 vibration modes were 
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used in the analysis shown in figure 3. The effect on both displacements and sensitivities 
of including fewer modes for each of the two load cases will be discussed below. 

Design Variables 

Two different classes of design variables were considered. The first design variable is a 
concentrated mass (initially zero) at the location of the translational damper. This design 
variable was also considered in [15]. The derivatives of the system mass and stiffness 
matrices with respect to this design variable are constant and zero, respectively. As a 
consequence, the derivatives of the system matrices required in the semi-analytical method 
can be calculated exactly by a simple forward difference operator. The beam thicknesses 
in each of the five spans were also design variables. Derivatives with respect to the five 
thickness design variables showed similar characteristics. Herein, results for derivatives 
with respect to the thickness in the rightmost span, h, along with derivatives with respect 
to the point mass, m, are presented. 

Effect of Finite Difference Step size on Derivatives 

The three methods described above for calculating sensitivities all rely on finite differ- 
ence operators at some stage in the algorithm. The forward and central difference methods 
rely on the operators in equations (8) and (9) to calculate derivatives of response quanti- 
ties. In the semi-analytical method the derivatives of the system matrices are calculated 
using the forward difference operator in equation (8). In all three cases the finite difference 
step size must be selected so that the operator provides a reasonable approximation to the 
derivative. If the step size is too large, the error due to truncating the series approximation 
of the derivative is large. If the step size is too small, the numerical condition error in 
performing the function evaluations (dynamic analyses) becomes large. 

To assess the effect of step size on the calculation of sensitivities for the five-span- 
beam, derivatives were calculated using the three methods with various step sizes. In this 
study the beam is undamped and the ramp loading is applied. All 26 vibration modes were 
included in the analysis. Figure 4a shows the derivative of displacement tt 2 at critical point 
number 5 with respect to the point mass design variable, m, as a function of step size. 
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As mentioned above, the derivatives of the system matrices with respect to this design 
variable can be calculated exactly in the semi-analytical method. As a result, the semi- 
analytical derivative is approximately constant for the six-order-of-magnitude change in 
step size shown in the figure. The central difference method uses the higher order operator 
and provides good accuracy over most of the step size range shown in the figure. The 
forward difference operator provides good accuracy with the smaller step sizes but begins 
to diverge earlier than the central difference method for the larger step sizes. 

Figure 4b shows the derivative of displacement ui at critical point number 5 with 
respect to the right-most span thickness, h, as a function of step size. In this case the 
system mass matrix is a linear function of this design variable and its derivative can be 
computed exactly. The system stiffness matrix is a cubic function of this design variable 
and its derivative can only be approximated by the forward difference operator. Still, the 
derivative computed by the semi-analytical method is very accurate except for the largest 
step size and is no worse for this case than the much more costly central difference method. 
Again, the forward difference operator results in substantial errors for the larger step sizes. 

Because this example has a relatively small number of degrees of freedom there is little 
condition error when small step sizes are used. To assess the effects of condition error which 
would occur for larger problems, the derivatives for the five-span-beam problem were also 
calculated on a DEC VAX computer. In the VAX calculations many of the floating point 
operations are performed using 32 bit precision compared with the 60 bit precision used in 
the studies described above. The derivative of displacement U 2 at critical point number 5 
with respect to the point mass is plotted as a function of finite difference step size in figure 
4c. Derivatives are calculated using the semi-analytical method, the central difference 
method, and the forward difference method. For the larger step sizes, the results from all 
four methods agree well with those calculated using 60 bit precision. For step sizes smaller 
than 10 -4 there is considerable error in the derivatives calculated using the forward and 
central difference methods. The derivatives calculated using the semi-analytical method, 
however, are highly accurate over the entire range of step sizes shown in the figure. 
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Modal Convergence of Displacements and Derivatives 

When vibration modes or other functions are used to reduce the basis in a transient 
response problem (equation 2), the key question is how many modes are required for an 
accurate solution. When the objective of the analysis is the calculation of sensitivities, 
small errors in the solution due to a reduced set of modes can have a pronounced effect on 
derivatives. 

Figure 5 shows the effect of number of vibration modes used in the analysis of the 
five-span-beam on selected displacements (shown by subscripts) at selected critical points 
(in parenthesis). Results for ramp and step loads acting on the undamped beam and 
the step load acting on the damped beam sure shown. The displacement /critical point 
combinations were selected to be representative of both the largest and smallest critical 
values. In all cases the convergence is very good with approximately ten modes yielding a 
converged solution. Although, as shown in figure 3a, there is considerable high frequency 
participation in the displacement history, the displacements at critical points are dominated 
by low frequency modes. 

The modal convergence of the displacement derivatives is considerably different, how- 
ever. Selected derivatives with respect to the point mass design variable are shown in 
figure 6 as a function of number of modes in the analysis. All derivative results shown 
in the figure were obtained using the central difference method with updated vibration 
modes to reduce errors associated with the finite difference operator and the modes. The 
consequences of keeping the vibration modes fixed instead of updating them are discussed 
below. For the step loading case with no damping, the derivatives of ui and «2 show 
especially poor convergence with more than 20 modes required. Small errors in the so- 
lutions for ux and U 2 due to modal truncation, which are unnoticeable in figure 5, are 
being magnified by the differentiation. For the ramp loading case with no damping, where 
the higher frequency modes play a less important role in the solution (see figure 3b) the 
convergence of the derivatives is nearly as good as for the displacements. 

For the case of the step loading with damping, the convergence of the derivatives is 
not as good as for the ramp loading case but better than for the case of step loading on 
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the undamped beam. The damping appears to improve convergence in the step loaded 
beam by filtering the high frequency modes from the displacement response. 

Figure 7 shows modal convergence results for selected displacement derivatives with 
respect to the thickness in the right-most span of the beam. The results are similar to those 
for the derivatives with respect to the point mass. The convergence for the derivatives in 
the step loaded, undamped beam is also especially slow. Approximately 24 of the total 
26 modes are required for converged derivatives. As before, the derivatives for the ramp 
loaded, undamped beam show good convergence. Convergence of the derivatives for the 
damped beam under the step load is also acceptable (i.e. convergence was achieved with 
about 10 modes). 

Fixed Versus Updated Modes in Derivative Calculations 

As mentioned above, the computational cost of updating the vibration modes for the 
perturbed analyses is substantial. The question of whether the modes from an initial 
design can be used in a finite-difference-based procedure to calculate sensitivities of the 
transient behavior has received considerable attention in the literature. In ref. [16] it 
was shown that there is a substantial difference in the derivatives of aircraft flutter speeds 
when fixed modes are used rather than the updated modes. In ref. [15], however, there 
was little difference in the derivatives of damping ratios for the five-span-beam when either 
fixed or updated modes were used. This was investigated here using the same five-span, 
undamped beam under the step load. As shown in figures 6 and 7 where the derivatives 
were calculated using the central difference method with updated modes, convergence with 
respect to the number of modes is very slow. Figure 8 shows the derivatives of selected 
displacements with respect to the right-most beam thickness calculated using forward 
difference procedures. Results with both fixed and updated vibration modes are shown. 
Again, the convergence as a function of number of modes is poor. However, for all three 
derivatives, the results are nearly the same for both the fixed mode and updated mode 
cases. 

A final study was performed to assess whether fixed modes can also produce accurate 
derivatives when the structure is substantially perturbed. For this case the point mass was 
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increased from zero to about 20% of the mass of the beam resulting in substantial frequency 
changes (see Table 1). The ramp loading was applied and the beam was undamped. The 
vibration mode shapes used in the forward difference method were calculated for the beam 
with zero point mass. Derivatives of selected beam displacements plotted as a function of 
the number of these modes are shown in figure 9. Convergence as a function of number of 
modes is fairly good in each case indicating that, at least for this example, the fixed modes 
can be used even after significant changes have been made to the original structure. 

Concluding Remarks 

Several computational concerns which arise in calculating the sensitivity of response 
quantities to design variable changes in structures under transient loading have been con- 
sidered. The structure can include general forms of damping and be acted on by general, 
transient loading. Natural vibration mode shapes are used to reduce the basis of the re- 
sulting transient problem. A critical point constraint formulation is used where explicit 
constraints are placed only on the largest values of a response quantity along a time history. 
Three sensitivity calculation methods are investigated in this study. Two of the methods 
involve application of the well-known one-point forward difference and two-point central 
difference operators. A third method involves explicit differentiation of the equations of 
motion to yield a set of differential equations in the sensitivities. The procedures are 
denoted forward difference, central difference, and semi-analytical respectively. All three 
methods use a highly accurate, explicit, time integration technique to reduce numerical 
errors in the integration of the equations of motion. 

The three methods are applied to a five-span damped beam acted on by an applied 
moment. Both a step and a ramp time distribution of the moment are considered. Two 
important issues investigated are the effect of finite difference step size on derivative ac- 
curacy and the effect of number of included vibration modes on derivative accuracy. The 
central difference method is relatively insensitive to finite difference step size but is com- 
putationally costly because two analyses of perturbed designs are required. Results from 
the forward difference method were shown to be accurate for small step sizes but in er- 
ror for larger step sizes. The semi-analytical method, on the other hand, was shown to 
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be as insensitive to step size as the central difference method but with a computational 
cost only slightly greater than the forward difference method. Calculations using all three 
methods performed on a computer with lower floating point precision indicate that the 
semi-analytical method may be more insensitive to condition errors than the forward or 
central difference methods. 

The convergence of derivatives from ail three methods as a function of number of 
modes in the analysis is similar. The time distribution of the loading, however, has a 
significant effect on modal convergence. Even though displacement quantities show good 
convergence for both the ramp and step loads, the convergence of derivatives for the beam 
under the step load is poor. For the five-span-beam example, the sensitivities calculated 
using fixed vibration modes are nearly identical to those calculated using updated modes. 
Using fixed modes substantially reduces the computational cost of sensitivity calculations. 
In addition, results from a computer with lower floating point precision indicate a potential 
improvement in the results for small step sizes when fixed rather than updated modes are 
used in the calculations. It was shown that damping improves the modal convergence. 

Based on insensitivity to finite difference step size, modal convergence, and compu- 
tational cost, the semi-analytical method appears to be the best choice for calculating 
sensitivities in transient response problems. 
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Figure 1.- Example displacement time history illustrating the critical point constraint selection 
process. 
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Figure 2.- Five-span-beam with applied end moment. 
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Figure 3.- Selec 
ment 
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Figure 3. Continued. 
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Figure 3.- Concluded. 



_Q 


_C 

CM 

D 

"D 


E 

■O 




CM 

'O 


CO 

'O 


'O 


LO 

'O 


CD 

'O 


'o 


28 


Finite difference stepsize, lb 

(a) Derivatives of displacement U 2 with respect to point mass 

Figure 4.- ElTect of finite difference stepsize on the accuracy of selected sensitivities (undamped, 
ramp loading). 




29 


(b) Derivative of displacement u\ with respect to rightmost span thickness. 
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Finite difference stepsize, in. 

(c) Derivative of displacement U 2 with respect to point mass calculated on the VAX 
computer. 

Figure 4.- Concluded. 



(2) , step load, undamped 

(3) , step load, undamped 
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Figure 5.- Modal convergence of selected displacements in the five-span-beam. 
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du 2 /dm (5), step load, undamped 
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Figure 6.- Modal convergence of displacement derivatives with respect to the point mass design 
variable for the five-span-beam (central difference method). 


du 4 /dh (5), step load, undamped 
du 1 /dh (5), step load, undamped 
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Number of modes 

Figure 7 .- Modal convergence of displacement derivatives with respect to thickness design vari- 
able for the five-span- be am (central difference method). 
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Figure 8.- Modal convergence of displacement derivatives calculated with both fixed and updated 
vibration modes (undamped, step load, forward difference methods). 


O du ^ /dm(3), ramp load 
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L igure 9. Modal convergence of displacement derivatives for the five-span-beam with a largi 
point mass (undamped, ramp load, forward difference method). 
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